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ABSTRACT:  Numerical  methods  based  on  unstructured  grids,  with  irregular  cells,  usually 
require  discrete  shape  functions  to  approximate  the  distribution  of  quantities  across  cells. 
For  control-volume  mixed  finite-element  methods,  vector  shape  functions  are  used  to  approx¬ 
imate  the  distribution  of  velocities  across  cells  and  vector  test  functions  are  used  to  minimize 
the  error  associated  with  the  numerical  approximation  scheme.  For  a  logically  cubic  mesh, 
the  lowest-order  shape  functions  are  chosen  in  a  natural  way  to  conserve  intercell  fluxes 
that  vary  linearly  in  logical  space.  Vector  test  functions,  while  somewhat  restricted  by  the 
mapping  into  the  logical  reference  cube,  admit  a  wider  class  of  possibilities.  Ideally,  an  error 
minimization  procedure  to  select  the  test  function  from  an  acceptable  class  of  candidates 
would  be  the  best  procedure.  Lacking  such  a  procedure,  we  first  investigate  the  effect  of 
possible  test  functions  on  the  pressure  distribution  over  the  control  volume;  specifically,  we 
look  for  test  functions  that  allow  for  the  elimination  of  intermediate  pressures  on  cell  faces. 
From  these  results,  we  select  three  forms  for  the  test  function  for  use  in  a  control- volume 
mixed  method  code  and  subject  them  to  an  error  analysis  for  different  forms  of  grid  irreg¬ 
ularity;  errors  are  reported  in  terms  of  the  discrete  L 2  norm  of  the  velocity  error.  Of  these 
three  forms,  one  appears  to  produce  optimal  results  for  most  forms  of  grid  irregularity. 


1  INTRODUCTION 

For  simulation  of  two-dimensional  (2-D)  flow  in  heterogeneous  porous  media,  it  has  been 
shown  that  mixed  methods,  and  in  particular  the  control- volume  mixed  finite-element 
(CVMFE)  methods,  are  often  the  most  efficient  methods  for  solving  for  the  velocity  field  (e.g., 
Durlofsky  (1994),  Cai  et  al.  (1997)).  In  this  report,  three  different  test  functions  are  evaluated 
for  use  in  the  the  CVMFE  method.  The  performance  of  the  test  functions  is  evaluated  using 
the  L2  norm  of  the  velocity  error  from  a  number  of  test  cases.  The  three-dimensional  (3-D) 
algorithm  described  herein  is  based  on  the  CVMFE  methodology  as  developed  by  Cai  et  al. 
(1997)  for  the  simulation  of  Darcian  flow  in  two  dimensions.  In  the  process  of  extending  the 
CVMFE  methodology  to  three  dimensions,  some  of  the  deficiencies  in  the  work  by  Cai  et  al. 
(1997)  were  investigated,  which  led  us  to  believe  that  a  different  test  function  might  improve 
the  performance  of  the  method.  In  addition,  it  has  been  discovered  that  the  vector  shape 
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discretization  with  vertex  locations  v0oo,  v100,  v0io,  v0oi,  Vn0,  Vi0i,  v0n  and  vm  indicated. 

functions  used  in  this  original  investigation  do  not  scale  up  to  three  dimensions  as  well  as 
we  had  hoped.  While  this  latter  problem  is  still  under  investigation,  it  is  noteworthy  as  it 
also  affects  the  interpretation  of  these  L2  results. 

In  the  CVMFE  method,  the  domain  is  discretized  into  potentially  irregular  hexahedral 
cells;  that  the  cells  can  have  irregular  shapes  allows  for  the  modeling  of  complex  hydrogeo¬ 
logical  systems.  Vector  shape  functions,  as  described  subsequently,  serve  as  basis  functions  to 
interpolate  the  velocity  over  the  interior  of  cells.  Vector  test  functions,  on  the  other  hand,  are 
used  as  weighting  factors  for  integrating  the  Darcy  relation  over  control  volumes  associated 
with  cell  faces;  this  usage  can  be  viewed  as  an  error  minimization  step  in  the  control  volume 
technique.  As  developed  for  CVMFE  methods,  an  additional  requirement  is  placed  on  the 
test  functions  in  that  their  use  is  expected  to  facilitate  the  elimination  of  pressures  on  faces 
between  cells;  these  pressures  are  extraneous  to  the  mixed-method  algorithm.  Our  intent  is 
to  develop  the  CVMFE  method  for  the  simulation  of  flow  in  3-D  heterogeneous  media. 

In  the  subsequent  three  sections  of  this  report,  background  information  concerning  the 
nature  of  the  problem  being  solved  and  the  solution  technique  are  discussed.  In  particular, 
in  section  3  the  vector  shape  functions  used  in  this  study  are  presented,  and  in  section  4 
the  rationale  for  the  three  test  functions  is  detailed.  The  first  test  problem,  for  uniform  flow 
in  a  simply  distorted  three-dimensional  domain,  follows  in  section  5.  In  section  6,  a  two- 
dimensional  nonuniform  flow  problem  is  presented  for  which  an  analytical  solution  exists.  A 
brief  discussion  of  the  results  concludes  the  report.  The  discussion  of  the  CVMFE  method 
presented  herein  is  necessarily  brief;  the  reader  is  invited  to  examine  the  work  by  Cai  et  al. 
(1997)  for  a  complete  explanation  of  the  2-D  algorithm. 

The  algorithm  used  in  this  study  is  strictly  three  dimensional;  when  2-D  simulations 
were  desired,  a  grid  was  constructed  from  a  single  layer  of  cells  of  uniform  thickness  such 
that  the  surface  of  this  layer  gives  the  desired  2-D  grid.  When  the  correct  test  function  is 
used,  this  procedure  results  in  matrix  equations  that  are  identical  to  those  resulting  from 
the  2-D  algorithm  of  Cai  et  al.  (1997).  Thus,  the  2-D  simulation  presented  herein  has  all  the 
characteristics  of  simulation  from  a  strictly  2-D  algorithm. 

2  BASIC  EQUATIONS 

The  numerical  method  outlined  in  this  paper  is  based  on  the  following  steady  flow  equation 
applied  within  a  3-D  domain  0: 

V-q  =  W(x,y,z)  (x,y,z)eQ  (1) 
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Here  q  is  the  specific  discharge  vector  and  W  (x,  y.  z )  is  a  source  term.  On  the  surface  Oil  of 
the  domain,  boundary  conditions  can  consist  of  fluxes  over  dQf  or  specified  pressures  over 
dilp.  The  specific  discharge  q  is  defined  from  the  Darcy  relation: 


q 


i> 


(2) 


where  p  is  the  pressure,  K(x,  y,  z )  is  the  hydraulic  conductivity  tensor  and  //  is  the  viscosity 
(the  gravitational  potential  is  neglected  from  (2)  for  notational  convenience).  For  the  CVMFE 
method,  the  hydraulic  conductivity  tensor  is  inverted,  causing  (2)  to  be  written  as 


Vp  =  -fiK  q 


(3) 


The  mixed-method  development  herein  uses  (1)  and  (3)  as  a  basis  of  the  numerical  approx¬ 
imation. 

The  domain  is  discretized  using  a  logically  regular  mesh  of  hexahedral  cells.  These  cells 
may  be  irregular  in  construction  to  the  point  where  faces  of  the  cells  are  not  planar.  Any 
given  hexahedral  cell  Q  can  be  described  from  the  location  of  its  vertices  at  v0oo5  Vioo,  v0io, 
vooi,  vno,  vioi,  von  and  vm,  where  vQ(g7  =  (xa3,n  yafjl,  za/3l)  (Figure  1).  A  trilinear  mapping 
exists  such  that  the  hexahedral  cell  Q  is  the  image  of  a  regular  reference  hexahedron,  Q, 
consisting  of  unit  cube  with  fixed  vertices  at  (0,0,0),  (1,0,0),  (0,1,0),  (0,0,1),  (1,1,0), 
(1, 0, 1),  (0, 1, 1)  and  (1, 1, 1)  on  a  point  by  point  basis;  Figure  1  depicts  these  relations.  This 
mapping  allows  that,  for  any  reference  location  r  =  (x,  y.  z)  within  Q,  the  equivalent  location 
r  =  (x,  y.  z)  within  Q  is  obtained  from  the  following  expression: 


r  =  v0  +  vax  +  vby  +  vcz  +  \dxy 


+vexz  +  Vfyz  +  Vgxyz 


(4) 


where  v0  -  v000,  va  -  vi00  -  v0,  v6  -  v0i0  -  v0,  vc  -  v00i  -  v0,  vd  -  vn0  -  v0  -  va  -  v6, 

ve  =  Vioi-V0-Va-Vc,  V;  =  V0n-Vo-V6-Vc,  Vg  =  Vm  -  V0  -  Va  -  Vb  -  Vc  -  Vd  -  Ve  -  V/. 

This  mapping  is  an  extension  of  the  conventional  bilinear  mapping  for  a  logically  rectangular 
grid  (Cai  et  al.  1997;  Garanzha  &;  Konshin  1999).  Note  that,  should  x  be  fixed  in  r,  then  a 
face  normal  to  x  within  Q  is  determined.  This  face  is  mapped,  via  (4),  into  an  equivalent  face 
within  Q.  Covariant  vectors,  defined  as  X(y,  z)  =  dr / dx ,  Y (x,  z)  =  dr/dy  and  Z(x,  y)  = 
dr /dz,  allow  for  definition  of  the  geometry  of  Q.  The  volumetric  Jacobian  J  for  passing  from 
Q  to  Q  is  simply  (Hildebrand  1962) 

J (x,  y,  z)  =  X(y,  z)  ■  (Y (x,  z)  x  Z(x,  y))  (5) 

while,  for  a  face  in  the  logical  x  direction,  the  surface  Jacobian  rx.  becomes  (Hildebrand 
1962) 

r X(x,y,z)  -  |Y(x,  z)  x  Z(x,  y)\  (6) 

Surface  Jacobians  can  similarly  be  defined  for  faces  in  the  logical  y  and  z  directions.  These 
few  relations  allow  us  to  define  the  necessary  quantities  used  within  the  CVMFE  method. 
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3  SHAPE  FUNCTIONS 

Velocity  basis  vectors  vt,j,k-m,n,p  are  associated  with  the  bulk  flux  through  each  face  of  a  cell 
Qi,j,k ;  here  %,  j,  k  represents  the  location  of  a  hexahedral  cell  within  the  discretized  domain. 
The  forms  of  the  vector  shape  functions  used  in  this  study  are 


V<j,fc;i/2,0,0  -  Pi+l/2,j,k  j 

Ji,j,k 


(7a) 


_  (1  J)^-i,j,k 

Vi,j,k;- 1/2, 0,0  -  Pi-l/2,j,k - j - 


Vij,k-,0, 1/2,0  —  Pi,j+ 1/2, k 
Vi,j,k;  0-1/2, 0  =  Pi,j-l/2,k 


0,1/2  —  Pi,j,  k+1/2 
0,0, -1/2  =  Pi,j,k- 1/2 


hj,k 


(1  —  zyZjjj; 


Jr 


hj,k 


(7b) 

(7c) 

(7d) 

(7e) 

(7f) 


where  is  the  volumetric  Jacobian,  and  are  the  covariant  vectors  for 

Qi,j,k,  and,  for  example,  pi+i/2,j,k  is  an  area  correction  factor  for  the  face  at  x  —  1  (denoted 
as  F,l+L/2.j,k)  such  that 


Pi+l/2,j,k{i)i  J) 


A+l/2, j,k 


Here,  ri+1/2j,k(t/,  5)  =  |Y  x  Z|i=1  is  the  surface  Jacobian  at  Fi+i/sj,k  and  Ai+ 1/2, j,k  is  the  area 
of  this  face.  Other  area  correction  factors  are  expressed  in  a  similar  manner.  These  shape 
functions  are  an  extension  of  those  used  by  Cai  et  al.  (1997)  for  simulation  of  flow  in  two 
dimensions.  Using  the  velocity  basis  vectors,  the  specific  discharge  can  be  approximated  at 
any  point  within  a  hexahedral  cell  Qt.3,k  as 


Q  ~  fi+l/2,j,k^i,j,k-,l/2,0,0  +  fi-l/2,j,kVi,j,k--l/2,0,0 


+fi,j+l/2,kVi,j,k-,0,l/2,0  +  fi,j-l/2,kVi,j-l/2,k-,0  -1/2,0 


+fiJ,k+l/2'Vhj,k-,0, 0,1/2  +  fi,j,k-l/2^i,j,k-,0,0,-l/2  (8) 

where,  for  example,  fi+i/2,j,k  is  the  bulk  flux  through  face  1/2, j,k-  If  this  approximation 
is  substituted  into  (1)  and  both  sides  are  integrated  over  Qi,j,ki  then  Gauss’s  divergence 
theorem  implies  that 


fi+l/2,j,k  fi-l/2,j,k  +  fi,j+l/2,k  fi.j-i/2.k 

+fi,j, k+1/2  -  fi,j,k- 1/2  =  /  w(x,  y,  z)  dxdydz  (9) 
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Figure  2:  Midsection  through  two  adjoining  cells  in  the  logical  x  direction.  Filled  and  open 
circles  represent  locations  on  faces  where  fluxes  are  estimated;  open  circles  represent  flux 
locations  on  faces  in  the  logical  y  direction  projected  onto  midsection.  The  x  symbol  is  used 
to  denote  pressure  locations  at  logical  cell  centers. 

Thus,  the  shape  functions  provide  for  a  discrete  form  of  the  continuity  equation  within 
the  discretized  region  0.  While  the  above  shape  functions  generally  appear  to  provide  an 
adequate  representation  of  the  specific  discharge  q  over  a  cell  in  two  dimensions,  their  use 
in  three  dimensions  is  more  questionable.  For  the  simulation  of  a  uniform  flow  field  in  three 
dimensions  with  an  irregular  grid,  it  can  be  shown  that  these  functions  do  not  always  properly 
capture  the  flux  distribution  across  every  irregular  cell.  This  result  may  explain  why,  even 
with  some  very  simple  irregular  grids,  exact  simulations  of  uniform  flow  are  not  achieved  in 
three  dimensions  with  the  present  CVMFE  algorithm. 

4  TEST  FUNCTIONS 

To  develop  a  discrete  form  for  the  Darcy  relation  (3),  a  control  volume  which  straddles  a  cell 
face  within  the  discretized  region  is  used;  a  control  volume  Q*+i/2,j,k  straddling  Fi+i/2,j,k  is 
depicted  in  Figure  2.  Weighted  with  a  test  function,  the  modified  Darcy  relation  (3),  with  an 
approximation  similar  to  (8)  for  the  specific  discharge,  is  integrated  over  a  straddle  control 
volume.  The  choice  of  the  test  function  is  motivated  by  the  efficiency  and  convenience  of 
eliminating  intermediate  pressures  associated  with  faces  adjoining  two  cells;  in  the  logical  x 
direction  of  Figure  2,  this  face  is  denoted  as  Fi+i/2,j,k-  A  simple  form  on  which  to  base  the 
test  function  in  the  logical  x  direction  is  the  covariant  vector  X;  one  form  suggested  by  Cai 
et  al.  (1997)  for  the  test  function  is 


where 


wi  • 


u 


on  Qi+ i/4j,* 
OH  Q*i+3/4,j,k 

otherwise 


Aj+i/4J,fc 


2 1  L  jl/2Ji^dxdydz 


Ai+3/4,j,* 


ml/2 

Ji+i,j,k  dxdydz 


(10) 
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This  test  function  shall  be  referred  to  herein  as  the  CJMR  function.  When  applied  to  the 
left-hand  side  of  (3),  this  test  function  gives  the  result 

/  Vp  •  wi+i/2,j,jfc  dxdydz 


1 

Aj+l/4  ,j,k 


f 

JQ 


i+l/4,j,k 


Vp  •  X^k  dxdydz 


+ 


A,-j 


.7. k  J  Q*: 


4+3/4J,k  JQ*+3/4,j,k 


Vp  •  Xj+1  dxdydz 


A, 


■*+1/4  ,j,fc 


r 

JQ* 


dp 


Qh  1/4J,*  ^ 


dxdydz 


+ 


A,; 


*+3/4,j,fc 


L 


dp 


q;+ 3/4J,* 


dxdydz 


1 

A*+l/4,j,fc 


+ 


1 


A,+3/4j,A 


dxdydz 


1 

— 

dp 

dp 
H-  -7— ' 

2 

z+i/4,j,/c 

z+3/4,j,/c_ 

—  Pi+l,j,k  Pi,j,k  (11) 

which  results  from  making  the  assumption  that  the  pressure  varies  linearly  in  reference  space 
Q  such  that 

dp/dx^+i/^jfi  =  2(pj_|_i//2J,A;  Pi,j,k)t 

(x,  y,  z)  e  Qi+i/4^ 

dp/dx  |*+3/4J,A  —  2(pj+ij^  —  Pi-\-l/2,j,k): 


(x,  y,  z)  €  Q-:3/,ij.a- 

For  an  irregular  grid,  even  in  the  presence  of  a  uniform  flow  field,  this  assumption  is  unlikely 
to  hold. 

An  alternate  test  function,  which  partially  repairs  this  problem,  is 


W?+l/2j',fc 


Xj5j.fc/A?:+i/4j^,  on  Qj+ i/4jyc 

^  '^*+ljjA/Aj-|-3/4J^:  Qi-\-S/4:,j,k 

^  0,  otherwise 


(12) 
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where  now 


w  =  L'L'L* 


A* 

1Vi+3/4,j,k 


1/2  i+1/4’j’k 
■  ^-i,j,kJi,j,k  dxdydz 

rl  ,1  y  1/2 

Jo  Jo  Jo  U?+3/4j’" 

■  -^*+i,j,fe'^(+i,j)fe  dxdydz 


and 


u 


*+1/4,  j,k 


u«+3/4,/,fc  ~ 


(  A^Cj-|-i/4,j,fc  j  ^l/i+l/4,/+’  ^+'i+l/4,/,fc) 

As2 

6*+l/4J,A 

(^i+3/4,j,fe,  ^^+3/4,/,*,  ^^+3/4,j,A:) 
As2 

^6i+3/4j,A 


and  As<+i/4j,*  =  (Ax2i+1/4Jjk  +  A yf+1/4dik  +  A ^+1/4J;fc)1/2,  Asi+3/4jJ-*  =  (Ax2i+3/^k 

+  AVi+3/4d,k  +  Azi+3/4,j,k)1/2-  Here:  (Axi+i/4,j,kj  AVi+i/4,j,k ,  A^+i/4,/,*)  is  the  vector  from  the 
physical  centroid  of  cell  Qyy  to  the  center  of  face  Fj+i/2jy;  (A®  i+3/4,j,k:  Alli+3/4,j,ki  Azi+3/4,j,k) 
goes  similarly  from  f++i/2jy  to  Qj+iyit-  When  applied  to  the  left-hand  side  of  (3),  this  test 
function  gives  the  result 


L  VP  ■  wi+i/2,j,k  dxdydz 

JQi+l/2,j,k 


-  -y—- - /  Vp  ■  X.(Jy  dxdydz 

^«+l/4j,A  J<^i+l/4,j,h 

+  t/— - -  /  Vp  •  x*+i,i,fc  dxdydz 

d\+3/4, j,k  JQ*+3/4 ,j,k 

= idr-  r  r  l  vp  ■  dMidi 


i+l/4,j,k 


d/2 

el  rl  rl/2 


1  r 1  r 1  /-V2 

+  Vp  •  Xyijy./yiyfc 

Ai+3/4,j,k  0  0  70 


~  Pi+l,j,k-  Pi,j,k 

The  last  result  follows  from  the  approximation  that 

dp  _  \  L.  i 

o/  \  ^  vP*,/,*  P*+l/2j,A;J  lu*+l/4,/,fch 

OP  u  j/  i  1/1J.A: 


(13) 
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(z,  y,  z)  G  Qi+i/4,j,k 


dp 

d(eu)i+3/4,j,k 


~  (pi+l/‘2,j,k  Pi+l,j,k )  |U'j+3/4J,fc|; 


(x,  y,  z)  G  Qi+3/4,j,k 

where  eu  is  the  unit  vector  u/|u|.  That  is,  the  pressure  variation  is  assumed  to  be  piecewise 
linear  in  physical  space  Q.  However,  this  procedure  does  not  insure  that  each  piecewise  linear 
approximation  in  the  logical  x.  y  and  z  directions  is  compatible  with  a  constant  Vp.  Because 
this  test  function  represents  a  modification  of  the  CJMR  test  function,  it  will  be  referred  to 
as  the  MCJMR  function. 

Recently,  Garanzha  and  Konshin  (1999)  have  proposed  a  simple  variant  of  the  test  func¬ 
tion;  their  suggested  form  is 


wi+l/2,j,k 


^i,j,k/  Ji,j,k,  On  1/4, j,k 

<  ^-i+l,j,k/  on  Qi+3/4,j,k 

k  0,  otherwise 


(14) 


This  test  function  shall  be  subsequently  referred  to  as  the  GK  function.  When  applied  to 
the  left  hand  side  of  (3),  this  function  yields 


/ 

Jq* 


Vp  •  w  **1/2^k  dxdydz 


=  [ 
Jq * 


^ P  ' 


Qi+l/4,j,k 


f 

Jq* 


®i+3/4,j,k  "A  I  l,j,k 

1  r  1  j- 1 
1/2 

1  rl  rl/2 


dxdydz 

dxdydz 


=  [  f  f  Vp  •  Xjj  *  dxdydz 

JO  Jo  J 1/2 

rl  rl  rl/2 

+  J  yo  yo  Vp-Xi+ij,*  dxdydz 


=  f  f  [  dxdydz 

Jo  Jo  j  1/2  OX  v/ 1  1/-1J.A- 

m1/2  <9p,  „  _  _ 

—  L*  dxdydz 

OX  Vi+3/4,l,fe 

=  /  l  (Pb=i-Pb=i/2)|«;+1/W<i»<<j 


9p. 


/  l  (/4  1/2  -  P|.r-o) 


+ 


+3/4,j,(! 


dydz 
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(15) 


~  P«+l ,j,k  ~  Pi,j,k 

Here,  the  major  approximation  is  that  the  integral  of  the  pressure  over  the  center  section 
in  reference  space  is  approximately  equivalent  to  the  average  pressure.  In  the  presence  of  a 
uniform  flow  field  and  in  two  dimensions,  this  approximation  is  exact  even  for  an  irregular 
grid;  in  three  dimensions,  even  in  the  presence  of  a  uniform  flow  field,  some  error  will  be 
introduced  when  an  irregular  grid  is  used  to  discretize  the  domain.  Overall,  this  test  func¬ 
tion  is  the  least  restrictive  as  far  as  assumptions  required  for  eliminating  the  pressure  on 
intermediate  faces. 

Over  the  straddling  control  volume  Q*+if2j,k  ^ie  specific  discharge  is  again  approximated 
as  a  linear  combination  of  fluxes  and  shape  functions  (7a)-(7/).  Allowing  q  to  be  this  ap¬ 
proximation,  then  at  any  point  in  Ql+1/2,j,k 

Q  =  fi+l/2,j,k  (vy,fc;  1/2, 0,0  +  Vj+ij;fe;_i/2,0,o) 

+fi-l/2,j,kVi,j,k--l/2,0,0  +  fi+3/2,j,kVi+l,j,k;l/2,0,0 

+fi,j+l/2,kVi,j,k-,0,l/2,0  +  fi,j-l/2,kVi,j,k-,0  -1/2,0 

+fi,j,k+l/2Vi,j,k-,0, 0,1/2  +  fi,j,k-l/2Vi,j,k-,0,0-l/2 

+fi+l,j+l/2,kVi+l,j,k$M2fi 

-1/2,0 

+fi+l,j,k+l/2'Vi+hi,k;0, 0,1/2 

+fi+l,j,k-l/2^i+l,j,k-,0, 0-1/2  (16) 

Assuming  that  any  of  the  above  forms  can  be  used  as  a  test  function,  then  a  discrete  form 
for  the  Darcy  relation  becomes 

Pi+i,j,k  ~  Pi,j,k  =  -fi  I  (K_1q)  •  w  dxdydz  (17) 

JQi+ k 

where  w  can  be  any  of  the  above  test  functions.  Completion  of  the  indicated  integrations  in 
(17)  results  in  a  set  of  coefficients  relating  the  bulk  fluxes  on  the  various  faces  of  Qi,j,k  and 
Qi+i,j,k  to  the  difference  in  pressures  at  the  cell  centers  (at  x  —  1/2,  y  —  1/2,  z  —  1/2  in 
each  cell).  This  method  is  applied  to  all  interior  faces  of  the  domain;  thus  a  discrete  Darcy 
relation  is  formed  for  every  interior  face.  These  relations  and  the  discrete  continuity  equation 
(9),  applied  to  every  cell  within  the  domain,  form  the  sets  of  equations  which  are  the  basis 


9 


Figure  3:  Method  by  which  distortion  is  introduced  into  an  otherwise  orthogonal  grid  con¬ 
sisting  of  4  x  4  x  4  cells;  vertex  at  center  of  domain  is  displaced  into  neighboring  cells  by 
distance  d. 

of  the  CVMFE  method.  In  applications  to  irregular  grids,  this  system  of  equations  is  non- 
symmetric;  otherwise,  it  is  similar  if  not  identical  in  structure  to  equations  that  result  from 
other  mixed-method  techniques.  The  Schur  complement  (Dougherty  1990)  was  used  to  solve 
these  equations  for  the  test  problems  presented  herein. 

5  THREE-DIMENSIONAL,  UNIFORM  FLOW  TESTS 

In  this  section,  the  effect  of  grid  distortion  on  the  discrete  L2  norm  of  the  velocity  error  is 
briefly  investigated.  Simulations  with  the  CVMFE  method  of  2-D  uniform  flow  on  irregular 
grids  give  the  result  that  the  fluxes  on  the  internal  faces  of  the  grid  are  estimated  essentially 
exactly  when  the  GK  test  function  is  used;  this  result  is  in  agreement  with  the  work  of 
Garanzha  and  Konshin  (1999).  As  there  is  always  some  error  associated  with  simulations 
using  the  CJMR  or  MCJMR  test  functions,  our  2-D  results  always  indicated  that  the  GK 
test  function  is  superior  for  simulating  uniform  flow  on  irregular  grids.  In  three  dimensions, 
simulations  of  uniform  flow  on  irregular  grids  with  the  GK  test  function  do  not  produce  exact 
flux  results  on  interior  faces.  We  attempt  to  quantify  the  difference  in  performance  between 
the  three  test  functions  using  a  simple  grid-distortion  test;  in  particular,  grid  distortion 
which  results  in  cells  with  non-planar  faces  is  considered.  A  test  module  was  created  by 
considering  uniform  flow  through  a  3-D  domain  consisting  of  4  x  4  x  4  cells  in  which  all  the 
cell  edges  are  equidimensional  and  orthogonal.  As  depicted  in  Figure  3,  a  distortion  in  the 
otherwise  regular  discretization  is  produced  by  moving  the  vertex  at  the  center  of  the  grid 
along  the  adjacent  cell  edge  boundary  in  the  x  direction;  the  four  faces  adjoining  this  vertex, 
which  were  originally  perpendicular  to  the  x  direction,  become  increasingly  non-planar  as  the 
offset  d  is  increased.  Constant  flux  boundary  conditions  are  imposed  at  both  the  upstream 
(x  =  0)  and  downstream  (x  =  4)  ends  of  the  domain  so  as  to  ensure  a  uniform  flow  field. 
The  numerical  model  produces  essentially  exact  flux  estimates  for  all  interior  faces  of  the 
discretized  domain  when  d  —  0;  that  is,  on  orthogonal  grids  all  three  test  functions  are 
equivalent  and  produce,  for  uniform  flow,  exact  simulations.  As  d  increases,  errors  in  the 
numerical  flux  estimates  are  to  be  found  on  all  faces  in  the  domain,  but  are  greatest  near 
the  distortion.  As  flow  is  uniform  and  in  the  x  direction,  the  normal  fluxes  across  all  faces 
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Figure  4:  L 2  norm  results  for  various  test  functions  from  the  uniform  flow,  simply  distorted 
grid  experiment.  Solid  line:  GK  test  function;  dashed  line:  MCJMR  test  function;  dashed 
and  dotted  line:  CJMR  test  function. 

are  easily  calculated  by  projecting  cell  facial  areas  onto  the  upstream  face  of  the  domain. 
The  exact  fluxes  were  compared  with  the  calculated  fluxes  using  L2  norms.  More  precisely, 
the  mean  normal  velocity  error  for  each  face  was  obtained  by  dividing  the  difference  of  these 
fluxes  by  the  facial  area;  then  a  root-mean-square  of  these  errors,  weighted  by  the  relative 
volume  of  the  control  volume  associated  with  each  face,  was  calculated.  This  is  essentially 
an  L 2  norm  of  the  velocity  error,  normalized  by  the  L2  norm  of  a  unit  vector. 

Figure  4  is  a  depiction  of  the  L 2  norm  results  as  a  function  of  the  relative  displacement 
d  for  the  three  test  functions  given  previously.  From  these  results,  the  performance  of  the 
algorithm  with  the  GK  test  function  is  shown  to  be  significantly  superior  to  that  of  the  CJMR 
or  MCJMR  test  functions;  the  CVMFE  method  with  the  GK  test  function  outperforms  the 
CJMR  test  function  by  an  order  of  magnitude.  Other  3-D  simulations  of  uniform  flow  in  the 
presence  of  an  irregular  grid  that  we  have  experimented  with  also  indicated  the  superiority 
of  the  GK  test  function;  however,  the  difference  in  performance  of  the  various  test  functions 
may  not  be  as  great.  For  example,  simulations  of  uniform  flow  have  been  performed  with 
grids  consisting  of  4  x  4  x  4  cells  in  which  all  the  internal  cell  faces  are  planar  but  have 
limited  random  orientations;  in  these  simulations,  the  exterior  faces  of  the  simulation  domain 
were  planar  and  orthogonal.  For  this  particular  example,  the  CVMFE  method  with  the 
GK,  MCJMR  and  CJMR  test  functions  give  L2  norms  of  1.04  x  10-3,  1.27  x  10-3  and 
1.46  x  10-3,  respectively.  However,  that  none  of  these  3-D  simulations  of  uniform  flow  on 
irregular  grids  with  random  faces  are  exact,  or  even  nearly  so,  causes  us  to  suspect  that  that 
the  aforementioned  problems  with  the  3-D  shape  functions  may  be  affecting  these  results.  If 
this  is  the  case,  then  these  L2  norm  results  may  largely  reflect  the  difficulty  of  using  shape 
functions  (7a)-(7/)  as  velocity  basis  functions  for  3-D  simulations.  Thus,  for  other  forms 
of  grid  irregularity  we  can  only  say,  at  this  time,  that  using  the  GK  test  function  in  the 
CVMFE  method  appears  to  be  marginally  superior  for  simulating  3-D  uniform  flow,  but  this 
marginal  performance  may  in  large  part  be  a  result  of  using  shape  functions  ( 7a)-(7f )  in  the 
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Figure  5:  Core  discretization  for  disk-shaped  heterogeneity  problem.  Shaded  area  represents 
region  of  conductivity  Ki;  remaining  domain  has  conductivity  K2.  Flow  is  largely  from  left 
to  right. 

CVMFE  algorithm.  An  example  of  an  essentially  exact  simulation  for  2-D,  uniform  flow  is 
given  in  the  next  section. 

6  TWO-DIMENSIONAL,  NONUNIFORM  FLOW  TEST 

In  order  to  examine  the  performance  of  the  various  test  functions  under  nonuniform  flow 
conditions,  a  2-D  problem  consisting  of  a  disk-shaped  heterogeneity  embedded  in  an  oth¬ 
erwise  uniform  medium  was  selected  (Figure  5).  In  this  test  problem,  nonuniform  flow  is 
produced  in  an  otherwise  uniform  flow  field  by  altering  the  hydraulic  conductivity  in  the 
embedded,  disk-shaped  heterogeneity  so  that  it  contrasts  with  the  hydraulic  conductivity  of 
the  remainder  of  the  domain.  Here,  the  hydraulic  conductivity  Ki  of  the  embedded,  disk¬ 
shaped  medium  is  assumed  to  be  equal  to  or  greater  than  the  hydraulic  conductivity  K2  of 
the  embedding  medium.  When  this  heterogeneous  geometry  is  posed  for  an  infinite  domain, 
an  analytical  solution  is  available  that  gives  pressure  as  a  function  of  space  and  the  contrast 
in  hydraulic  conductivity  between  the  disk-shaped  heterogeneity  and  the  remainder  of  the 
domain  (Carslaw  &  Jaeger  1959);  in  this  case,  flow  is  assumed  to  be  uniform  at  an  infinite 
distance  from  the  heterogeneity.  As  depicted  in  Figure  5,  this  heterogeneous  geometry  is 
relatively  easy  to  discretize  for  a  finite  domain;  here,  it  is  assumed  that,  in  the  absence  of 
a  contrast  in  hydraulic  conductivity,  flow  would  be  uniform  and  from  left  to  right  in  the 
x  direction.  For  simulation  purposes,  the  core  discretization  depicted  in  this  figure  was  ex¬ 
tended  by  adding  four  more  layers  of  cells  at  the  x  =  0,  x  =  8  and  z  =  4  boundaries;  these 
cells  were  all  orthogonal  with  the  dimensions  1/2  x  1/2.  The  purpose  of  these  additional 
cell  layers  was  to  shift  the  boundary  conditions  farther  from  the  disk-shaped  heterogeneity 
and  thus  minimize  their  influence  on  the  resulting  simulations.  Flux  boundary  conditions 
were  imposed  at  x  =  —2,  x  =  10  and  z  =  6;  these  boundary  conditions  were  extracted 
from  the  analytical  solution  by  integrating  the  normal  velocity  component,  as  obtained  from 
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Figure  6:  L 2  norm  results  for  various  test  functions  from  the  2-D,  nonuniform  flow  experiment. 
Solid  line:  GK  test  function;  dashed  line:  MCJMR  test  function;  dotted  line:  CJMR  test 
function. 

the  analytical  solution,  over  the  exterior  cell  faces  on  these  boundaries.  For  each  different 
hydraulic  conductivity  contrast,  it  was  necessary  to  generate  a  new  flux  boundary  condition. 
Facially  averaged  velocities  from  the  analytical  solution  were  used  as  exact  velocities  in  cal¬ 
culating  the  L 2  norms;  the  L2  norms  were  calculated  only  over  the  irregular  faces  of  the  core 
discretization  shown  in  Figure  5. 

The  performance  of  the  various  test  functions  for  this  test  case,  as  described  by  the  L2 
norm,  is  depicted  in  Figure  6.  Here,  the  L 2  norm  results  are  presented  as  a  function  of  the 
ratio  Ki/K2,  representing  the  contrast  in  hydraulic  conductivity,  which  was  allowed  to  vary 
from  unity  (no  contrast)  to  ten  (high  contrast).  When  Ki/K2  =  1,  the  L2  norm  for  the 
simulations  using  the  GK  test  function  is  null,  while  the  simulations  with  the  CJMR  or 
MCJMR  test  functions  are  finite  and  moderately  significant;  this  then  is  an  example  of  2-D 
uniform  flow  simulation  where  the  results  from  using  the  GK  test  function  are  exact.  However, 
the  difference  in  L2  norms  for  simulations  with  these  various  test  functions  decreases  rather 
rapidly  as  the  ratio  K\JK2  increases;  when  Ki/K2  =  10,  there  is  little  difference  between  the 
various  L 2  norm  results.  Apparently,  as  nonuniform  flow  begins  to  dominate  the  simulation, 
it  becomes  more  difficult  for  the  shape  functions  to  locally  approximate  the  velocity  field  and 
the  simulation  error  is  increased.  Thus,  when  flow  is  strongly  nonuniform,  the  approximation 
associated  with  the  shape  functions  (7a)-(7f)  appears  to  be  a  greater  source  of  error  than 
that  associated  with  the  test  function. 

7  CONCLUDING  REMARKS 

In  general,  results  given  here  indicate  that  the  GK  test  function  proposed  by  Garanzha  and 
Konshin  (1999)  produces  the  best  overall  results.  While  it  is  possible  to  simulate  exactly 
uniform  flow  on  two-dimensional,  irregular  grids  using  the  GK  test  function,  we  have  not 
been  able  to  repeat  this  feat  in  three  dimensions.  This  may  be  the  result  of  a  poor  scaling 
of  the  shape  functions  between  two  and  three  dimensions,  but  may  also  result,  in  part,  from 
a  difference  in  the  averaging  properties  of  the  test  function.  That  is,  while  the  GK  test 
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function  in  two  dimensions  gives  a  true  average  of  the  pressure,  the  average  is  weighted  in 
three  dimensions.  We  suspect,  however,  that  deficiencies  in  the  shape  functions  ( 7a)-(7f ) 
are  a  more  significant  factor  in  the  inability  of  the  current  CVMFE  algorithm  to  exactly 
simulate  3-D  uniform  flow  on  irregular  grids.  In  the  simulation  of  a  highly  nonuniform,  two- 
dimensional  flow  field,  the  superiority  of  the  GK  test  function  is  less  evident;  this  may  simply 
be  the  result  of  using  a  rather  coarse  grid  for  the  simulation  in  this  test.  The  GK  test  function 
does  appear  to  perform  better  in  simulations  involving  non-planar  faces  on  cells.  Non-planar 
faces  are  expected  to  be  a  common  feature  of  3-D,  irregular  grids;  thus,  the  GK  test  function 
is  likely  to  be  an  important  factor  in  extending  the  CVMFE  methology  to  three  dimensions. 
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